# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt

X = np.empty( ( 100, 2 ) )
X[:, 0] = np.random.uniform( 0, 100, size = 100 )
X[:, 1] = 0.75 * X[:, 0] + 3. + np.random.normal( 0, 10, size = 100 ) 

plt.scatter( X[:, 0], X[:, 1] )
plt.show()

def demean( X ):
    return X - np.mean( X, axis = 0 )
X_demean = demean( X )
plt.scatter( X_demean[:, 0], X_demean[:, 1] )
plt.show()

# 公式
def f( w, X ):
    return np.sum( X.dot( w ) ** 2 ) / len( X )

# 偏导数
def df( w, X ):
    return X.T.dot( X.dot( w ) ) * 2 / len( X )

def direction( w ):
    return w / np.linalg.norm( w )  #该函数可求出向量的模

# 计算w的方向
def first_component( df, X, initial_w, eta, n_iters = 1e4, epsilon = 1e-8 ):
    w = direction( initial_w ) # 初始化进来的w，进行向量化
    cur_iter = 0
    
    while cur_iter < n_iters:
        gradient = df( w, X )
        last_w = w
        w = w + eta * gradient
        w = direction( w )
        
        if ( abs( f( w, X ) - f( last_w, X ) ) ) < epsilon:
            break
        
        cur_iter += 1
    
    return w
    
initial_w = np.random.random( X.shape[1] )
eta = 0.01
w = first_component( df, X, initial_w, eta )
print( w )  # 第一主成分的方向确定

plt.scatter( X_demean[ :, 0 ], X_demean[ :, 1 ] )
plt.plot( ( 0, w[0]*50 ), ( 0, w[1]*50 ), color = "r" )
plt.show()

# 求第二主成分
''' 可以这样写，但是可以使用向量化
X2 = np.empty( X.shape )
for i in range( len( X ) ):
    X2[i] = X[i] - X[i].dot( w ) * w
'''
X2 = X - X.dot( w ).reshape( -1, 1 ) * w
    
plt.scatter( X2[ :, 0 ], X2[ :, 1 ] )
plt.show()

w2 = first_component( df, X2, initial_w, eta )
print( w2 )

print( w.dot( w2 ) ) # 约等于 0
    
# 封装一个函数，实现获得前n个主成分
def first_n_components( n, X, eta = 0.01, n_iters = 1e4, epsilon = 1e-8 ):
    X_pca = X.copy()
    X_pca = demean( X_pca )
    res = []
    for i in range( n ):
        initial_w = np.random.random( X_pca.shape[1] )
        w = first_component( df, X_pca, initial_w, eta )
        res.append( w )
        
        X_pca = X_pca - X_pca.dot( w ).reshape( -1, 1 ) * w
        
    return res

print( first_n_components( 2, X ) ) # 求出了前两个主成分 
#[array([ 0.77338545,  0.63393608]), array([-0.6339305 ,  0.77339002])]